home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
ada
/
gwuada_9.zip
/
ARITH.C
< prev
next >
Wrap
C/C++ Source or Header
|
1993-07-27
|
47KB
|
2,081 lines
/*
* Copyright (C) 1985-1992 New York University
*
* This file is part of the Ada/Ed-C system. See the Ada/Ed README file for
* warranty (none) and distribution info and also the GNU General Public
* License for more details.
*/
/* Avoid problems with library routines on PC that require pointer
* to double instead of double (this is violation of ANSI specs).
* for now do not compile these modules
*/
/*
* 11-oct-85 shields
* fix so can handle most negative integer properly. The changes
* are more in the nature of a patch than a solution. The proper
* way is NOT to convert negative values to positive form before
* processing; instead, positive values should be converted to
* negative and conversions done on negative values. However, this
* can be put off until later.
*
* 5-jun-85 shields
* add rat_tol and int_tol which are like rat_toi and int_toi, resp.,
* except that they return long results. These are needed to support
* fixed-point in generator.
*
* 1-aug-85 shields
* add calls to int_free() to free intermediate values known to be dead.
* add some copies where needed in building new rational values.
*/
#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include "config.h"
#include "miscp.h"
typedef struct Rational_s {
int *rnum; /* numerator */
int *rden; /* denominator */
} Rational_s;
typedef Rational_s * Rational;
#include "arithp.h"
static int *int_frl(long);
static int *int_ptn(int);
static int *subu_int(int *, int *);
static void int_div(int *, int *, int **, int **);
static int *int__addu(int *, int *);
static int *int__norm(int *);
static int *int_new(int);
static void int_free(int *);
static void rat_free(Rational);
static double pow2(int);
/* Some of the procedures want to signal overflow by returning the
* string 'OVERFLOW'. In C we do this by setting global arith_overflow
* to non-zero if overflow occurs, zero otherwise
*/
int arith_overflow = 0;
int ADA_MIN_INTEGER;
int ADA_MAX_INTEGER;
int *ADA_MAX_INTEGER_MP;
int *ADA_MIN_INTEGER_MP;
long ADA_MIN_FIXED, ADA_MAX_FIXED;
int *ADA_MIN_FIXED_MP, *ADA_MAX_FIXED_MP;
/* _LONG form is form that can be held in C (long) */
#ifdef MAX_INTEGER_LONG
long MIN_INTEGER_LONG;
int *MAX_INTEGER_LONG_MP;
int *MIN_INTEGER_LONG_MP;
#endif
#define ABS(x) ((x)<0 ? -(x) : (x))
#define SIGN(x) ((x)<0 ? -1 : (x) == 0 ? 0 : 1)
/*
* M U L T I P L E P R E C I S I O N I N T E G E R
*
* A R I T H M E T I C P A C K A G E
*
* Robert B. K. Dewar
* June 16th, 1980
*
* This package of routines implements multiple precision integer
* arithmetic using what are called the "classical algorithms" in
* Knuth "Art of Programming", Volume 2, Section 4.3.1. The style
* of the algorithms follows Knuth fairly closely, and section
* 4.3.1 can be consulted for further theoretical details.
*
* Multiple precision integers are represented as tuples of small
* integers in the range 0 to BAS - 1, where BAS is a power of 10
*(actually 10 ** DIGS) which is the base used to represent the
* long integers. Essentially the successive elements of the tuple
* are the digits of the representation in base BAS. All integers
* are normalized so that the first digit is non-zero, except in
* the case of zero itself. The sign is carried with the first
* digit value, all remaining digits are always positive.
*
* Some examples of the representation, assuming DIGS = 4 and
* BAS = 10000(note that the choice of BAS as a power of 10
* is for convenience of the conversion routines, and is not
* required for correct operation of the arithmetic algorithms).
* -123456789 [-1 , 2345 , 6789]
* 0 [0]
* 123456789 [1 , 2345, 6789]
* The constants BAS and DIGS must be defined as global constants
* in a program using these routines. It is assumed that the value
* BAS * BAS - 1 can be represented as a SETL integer value.
* The following routines are provided:
* int_abs integer absolute value
* int_add integer addition
* int_div integer division
* int_eql integer equal to
* int_exp raise multiple integer to multiple integer power
* int_fri convert multiple precision integer from integer
* int_frs convert multiple precision integer from string
* int_geq integer greater than or equal to
* int_gtr integer greater than
* int_len number of digits in multiple precision integer
* int_leq integer less than or equal to
* int_lss integer less than
* int_mod integer modulus
* int_mul integer multiplication
* int_neq integer not equal to
* int_ptn integer power of ten
* int_quo integer quotient
* int_rem integer remainder
* int_sub integer subtraction
* int_toi convert integer to C integer
* int_tos integer to string
* int_umin integer unary minus
*/
/* Internal procedures have names starting with int__ */
int *int_abs(int *u) /*;int_abs*/
{
/* Absolute value of multiple precision integer */
int *pu;
pu = int_copy(u);
if (pu[1] < 0)
pu[1] = -pu[1];
return pu;
}
int *int_add(int *u, int *v) /*;int_add*/
{
/* Add signed integers
* This procedure implements the familiar algorithm of comparing
* the signs, if the signs are the same, then the result is the
* sum of the magnitudes with the sign of the operands. If the
* signs differ, then the number with the smaller magnitude is
* subtracted from the larger magnitude and the result sign is
* that of the operand with the larger magnitude.
*/
int *t;
if (u[1] >= 0 && v[1] >= 0)
return (int__addu(u, v));
else
if (u[1] < 0 && v[1] < 0) {
u[1] = -u[1];
v[1] = -v[1];
t = int__addu(u, v);
u[1] = -u[1];
v[1] = -v[1];
t[1] = -t[1];
return t;
}
else {
int us, vs;
if (us = (u[1] < 0)) {
u[1] = -u[1];
}
if (vs = (v[1] < 0)) {
v[1] = -v[1];
}
if (int_gtr(u, v)) {
t = subu_int(u, v);
if (vs) v[1] = -v[1];
if (us) {
u[1] = -u[1];
t[1] = -t[1];
}
return t;
}
else {
t = subu_int(v, u);
if (us) u[1] = -u[1];
if (vs) {
v[1] = -v[1];
t[1] = -t[1];
}
return t;
}
}
}
int int_eql(int *u, int *v) /*;int_eql*/
{
/* Compare multiple precision integers for equality */
int n;
if (u[0] != v[0])
return FALSE;
for (n = u[0]; n > 0; n--)
if (u[n] != v[n]) return FALSE;
return TRUE;
}
int *int_exp(int *u, int *v) /*;int_exp*/
{
/* Raise a multiple precision integer to a multiple precision integer
* power using a modified version of the Russian Peasant algorithm
* for exponentiation.
*/
int sn;
int u1;
int i;
int carry;
int half;
int *running, *runningp;
int *result, *resultp;
/* Compute sign of result: positive if v is even, the sign of u if
* v is odd.
*/
sn =((v[v[0]] % 2) == 1) ? SIGN(u[1]) : 1;
u1 = u[1];
if (u[1] < 0)
u[1] = -u1;
v = int_copy(v);
/* Starting the result at 1 and running at u, loop through the binary
* digits of v, squaring running each time, and multiplying the result
* by the current value of running each time a 1-bit is found.
*/
result = int_con(1);
running = int_copy(u);
while(int_nez(v)) {
/* If v is odd then multiply result by running. */
if (v[v[0]] % 2 == 1) {
resultp = result; /* save current value */
result = int_mul(result, running);
int_free(resultp); /* free dead value */
}
/* Square running. */
runningp = running; /* save current value */
running = int_mul(running, running);
int_free(runningp); /* free dead value */
/* Halve v. */
carry = 0;
for (i = 1; i <= v[0]; i++) {
half = BAS * carry + v[i];
carry = half % 2;
v[i] = half / 2;
}
v = int__norm(v);
}
int_free(running);
int_free(v);
if (sn < 1)
result[1] = -result[1];
u[1] = u1;
return result;
}
int *int_fri(int i) /*;int_fri*/
{
/* Convert an integer to a multiple precision integer.
*
* First check the sign of i.
*/
int sn = 0;
int n = i;
int *t;
int ti;
int d;
if (i < 0) {
/* Special test for most negative as code below won't work